/* (c) Copyright 1993 by Joseph Felsenstein.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#include <boost/shared_ptr.hpp>
#include <boost/shared_array.hpp>
#include <math.h>
#include <iostream>
#include <sstream>
#include "KITSCH_nearest_neighbor.h"
#include "tree_to_newhampshire_factory.h"
#include "tree_to_paralogy_factory.h"
using std::vector;
using std::string;
using std::cerr;
using std::cout;
using std::ostringstream;
using boost::shared_array;
using paralogy_tree::t_tree;
using paralogy_tree::t_node;

using boost::shared_ptr;
using boost::shared_array;

namespace{
struct t_debug_function
{
	string					func_name;
	unsigned				priority;
	static std::ostream&	ostrm;
	static unsigned 		recurse;
	static unsigned 		indent;
	static unsigned 		debug_level;
	t_debug_function(unsigned i, const string& s): func_name(s), priority(i)
	{
		if (priority <= debug_level)
			ostrm << string(indent * recurse, ' ') << func_name << " entry\n";
		++recurse;
	}
	~t_debug_function()
	{
		--recurse;
		if (priority <= debug_level)
			ostrm << string(indent * recurse, ' ') << func_name << " exit\n";
	}
};

unsigned t_debug_function::recurse = 0;
unsigned t_debug_function::indent = 4;
unsigned t_debug_function::debug_level = 0;
std::ostream&	t_debug_function::ostrm = cerr;
#define debug_func(ii) t_debug_function debug_function(ii, __func__)
}

namespace kitsch_tree
{
//	for (unsigned i = 1; i < cnt_leaves; ++i)
//		for (unsigned j = 0; j < i; ++j)
//	for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
//		for (unsigned j = i + 1; j < cnt_leaves; ++j)
//	two forms are identical
	
	




#define epsilon         0.000001   /* a very small but not too small number */


struct kitsch_node
{		  /* describes a tip species or an ancestor */
	kitsch_node(unsigned count) : next (0), back(0), index(0), tip(false), d(count, 0.0),
						w(count, 1.0), t(0.0), sametime(false), weight(0.0), processed(false)
	{
	}
	kitsch_node *next;
	kitsch_node *back;
	unsigned index;
	bool tip;				/* present species are tips of tree       */
	vector<double> d;
	vector<double> w;		/* distances and weights                  */
	double t;				/* time                                   */
	bool sametime;			/* bookkeeps scrunched nodes              */
	double weight;			/* weight of kitsch_node used by scrunch         */
	bool processed;			/* used by evaluate                       */
	kitsch_node& operator = (const kitsch_node& other)
	{
		d			= other.d;
		w			= other.w;
		t			= other.t;
		sametime	= other.sametime;
		weight		= other.weight;
		processed	= other.processed;
		return *this;
	}
	

};

}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// globals
namespace
{

	vector<shared_ptr<kitsch_tree::kitsch_node> >	shared_mem;
	vector<kitsch_tree::kitsch_node* > 				treenode, bestree;
	vector<unsigned>									enterorder;
	kitsch_tree::kitsch_node*						root		= 0;
	kitsch_tree::kitsch_node*						best_root	= 0;
	// cnt_seq = number of species
	unsigned										cnt_seq;
	unsigned										cnt_nodes;
	unsigned	jumb		= 0;
	
	unsigned	njumble				= 10;
	double		power				= 2.0;
	unsigned	cnt_trees_examined	= 0;
	double		like				= 0.0;
	double		bestyet				= 0.0;
	kitsch_tree::kitsch_node *	there = 0;
}





namespace kitsch_tree
{

//________________________________________________________________________________________

//	doinit

//________________________________________________________________________________________
void doinit()
{
	/* initializes variables */

	treenode.insert(treenode.begin(), cnt_nodes, 0);
	for (unsigned i = 0; i < cnt_seq; ++i)
	{
		shared_mem.push_back(shared_ptr<kitsch_node>(new kitsch_node(cnt_nodes)));
		treenode[i] = shared_mem.back().get();
	}
	for (unsigned i = cnt_seq; i < cnt_nodes; ++i)
	{
		kitsch_node *q = 0;
		kitsch_node *p = 0;
		for (unsigned j = 0; j < 3; ++j)
		{
			shared_mem.push_back(shared_ptr<kitsch_node>(new kitsch_node(cnt_nodes)));
			p = shared_mem.back().get();
			p->next = q;
			q = p;
		}
		p->next->next->next = p;
		treenode[i] = p;
	}
	//cerr << "init 1\n";
	if (njumble > 1)
	{
		bestree.insert(bestree.begin(), cnt_nodes, 0);
		for (unsigned i = 0; i < cnt_seq; ++i)
		{
			shared_mem.push_back(shared_ptr<kitsch_node>(new kitsch_node(cnt_nodes)));
			bestree[i] = shared_mem.back().get();
		}
		for (unsigned i = cnt_seq; i < cnt_nodes; ++i)
		{
			kitsch_node *q = 0;
			kitsch_node *p = 0;
			for (unsigned j = 0; j < 3; j++)
			{
				shared_mem.push_back(shared_ptr<kitsch_node>(new kitsch_node(cnt_nodes)));
				p = shared_mem.back().get();
				p->next = q;
				q = p;
			}
			p->next->next->next = p;
			bestree[i] = p;
		}
		//cerr << "init 2\n";
	}
	//cerr << "init finished\n";
}  /* doinit */


void assigndata(std::vector<boost::shared_array<double> >& matrix,
					unsigned cnt_seq, unsigned cnt_nodes)
{
	debug_func(1);
	
	for (unsigned i = 0; i < cnt_nodes; ++i)
		treenode[i]->index = i;
	
	for (unsigned i = 0; i < cnt_seq; ++i)
		treenode[i]->tip = true;
	
	for (unsigned i = cnt_seq; i < cnt_nodes; ++i)
	{
		kitsch_node* p = treenode[i]->next;
		while (p != treenode[i])
		{
			p->index = i;
			p = p->next;
		}
	}
	
	
	//cerr << "assigndata1\n";
	
	
	for (unsigned i = 0; i < cnt_nodes; ++i)
		bestree[i]->index = i;

	for (unsigned i = 0; i < cnt_seq; ++i)
		bestree[i]->tip = true;

	for (unsigned i = cnt_seq; i < cnt_nodes; ++i)
	{
		kitsch_node* p = bestree[i]->next;
		while (p != bestree[i])
		{
			p->index = i;
			p = p->next;
		}
	}
	//cerr << "assigndata2\n";
		
		
	for (unsigned i = 0; i < cnt_seq - 1 ; ++i)
		for (unsigned j = i + 1; j < cnt_seq; ++j)
		{
			if (matrix[i][j] < epsilon)
				treenode[i]->d[j] = 
				treenode[j]->d[i] = epsilon;
			else
				treenode[i]->d[j] = 
				treenode[j]->d[i] = matrix[i][j];
			if (matrix[i][j] > 5000.0 || matrix[i][j] < 0.0)
				treenode[i]->w[j] = 
				treenode[j]->w[i] = 0.0;
			else
				treenode[i]->w[j] = 
				treenode[j]->w[i] =	1.0 / exp(2.0 * log(treenode[i]->d[j]));
		}
	
	
	
}




//________________________________________________________________________________________

//	add

//________________________________________________________________________________________
void add(kitsch_node* below, kitsch_node* newtip, kitsch_node* newfork)
{
	debug_func(1);
	
	/* inserts the nodes newfork and its left descendant, newtip,
	   to the tree.  below becomes newfork's right descendant */
	assert(below);
	assert(newtip);
	assert(newfork);
	if (below != treenode[below->index])
		below = treenode[below->index];
	//cerr << "add1\n";
	if (below->back != 0)
		below->back->back = newfork;
	newfork->back = below->back;
	below->back = newfork->next->next;
	newfork->next->next->back = below;
	newfork->next->back = newtip;
	newtip->back = newfork->next;
	//cerr << "add2\n";
	if (root == below)
		root = newfork;
	root->back = 0;
}  /* add */




//________________________________________________________________________________________

//	re_move

//________________________________________________________________________________________
void re_move(kitsch_node *& item, kitsch_node*& fork)
{
	debug_func(1);
	
	/* removes nodes item and its ancestor, fork, from the tree.
	   the new descendant of fork's ancestor is made to be
	   fork's second descendant (other than item).  Also
	   returns pointers to the deleted nodes, item and fork */

	if (item->back == 0)
	{
		fork = 0;
		return;
	}
	fork = treenode[item->back->index];
	if (root == fork)
	{
		if (item == fork->next->back)
			root = fork->next->next->back;
		else
			root = fork->next->back;
	}
	kitsch_node *p = item->back->next->back;
	kitsch_node *q = item->back->next->next->back;
	if (p != 0)
		p->back = q;
	if (q != 0)
		q->back = p;
	fork->back = 0;
	p = fork->next;
	while (p != fork)
	{
		p->back = 0;
		p = p->next;
	}
	item->back = 0;
}  /* remove */




//________________________________________________________________________________________

//	scrunchtraverse

//________________________________________________________________________________________
void scrunchtraverse(kitsch_node *u, kitsch_node** closest, double *tmax)
{
	//cout << "scrunch traverse tmax=" << *tmax << "\n";
	/* traverse to find closest kitsch_node to the current one */
	if (!u->sametime)
	{
		if (u->t > *tmax)
		{
			*closest = u;
			*tmax = u->t;
		}
		return;
	}
	u->t = treenode[u->back->index]->t;
	if (!u->tip)
	{
		scrunchtraverse(u->next->back, closest, tmax);
		scrunchtraverse(u->next->next->back, closest, tmax);
	}
}  /* scrunchtraverse */




//________________________________________________________________________________________

//	combine

//________________________________________________________________________________________
void combine(kitsch_node& a, kitsch_node& b)
{
	debug_func(4);
	/* put kitsch_node b into the set having the same time as a */
	if (a.weight + b.weight <= 0.0)
		a.t = 0.0;
	else
		a.t = (a.t * a.weight + b.t * b.weight) / (a.weight + b.weight);
	a.weight += b.weight;
	b.sametime = true;
}  /* combine */




//________________________________________________________________________________________

//	scrunch

//________________________________________________________________________________________
void scrunch(kitsch_node& s)
{
	debug_func(3);
	/* see if nodes can be combined to prevent negative lengths */
	kitsch_node *closest	= 0;
	do
	{
		double tmax	= -1.0;
		if (!s.tip)
		{
			scrunchtraverse(s.next->back, &closest, &tmax);
			scrunchtraverse(s.next->next->back, &closest, &tmax);
		}
		// found = (tmax > s.t);
		if (tmax <= s.t)
			break;
		combine(s, *closest);
	} while (1);
}  /* scrunch */




//________________________________________________________________________________________

//	secondtraverse

//________________________________________________________________________________________
void secondtraverse(kitsch_node* a, 
					kitsch_node* q, 
					kitsch_node* u, 
					kitsch_node* v, 
					unsigned i, unsigned j, unsigned k, double *sum)
{
	debug_func(4);
	
	/* recalculate distances, add to sum */

	if (!(a->processed || a->tip))
	{
		secondtraverse(a->next->back, q,u,v,i,j,k,sum);
		secondtraverse(a->next->next->back, q,u,v,i,j,k,sum);
		return;
	}
	if (!(a != q && a->processed))
		return;
	unsigned l = a->index;
	double wil = u->w[l];
	double wjl = v->w[l];
	double wkl = wil + wjl;
	double wli = a->w[i];
	double wlj = a->w[j];
	double wlk = wli + wlj;
	q->w[l] = wkl;
	a->w[k] = wlk;
	if (wkl <= 0.0)
		q->d[l] = 0.0;
	else
		q->d[l]	= (wil * u->d[l] + wjl * v->d[l]) / wkl;
	if (wlk <= 0.0)
		a->d[k] = 0.0;
	else
		a->d[k]	= (wli * a->d[i] + wlj * a->d[j]) / wlk;
	if (wkl > 0.0)
	{
		double TEMP = u->d[l] - v->d[l];
		(*sum) += wil * wjl / wkl * (TEMP * TEMP);
	}
	if (wlk > 0.0)
	{
		double TEMP = a->d[i] - a->d[j];
		(*sum) += wli * wlj / wlk * (TEMP * TEMP);
	}
}  /* secondtraverse */




//________________________________________________________________________________________

//	firstraverse

//________________________________________________________________________________________
void firstraverse(kitsch_node *q_, kitsch_node* r, double *sum)
{
	debug_func(3);
  /* firsttraverse                              */
   /* go through tree calculating branch lengths */
	/* Local variables for firstraverse: */
	//cout << "f sum=" << *sum << "\n";
	if (q_ == 0)
		return;
	
	kitsch_node *q = q_;
	q->sametime = false;
	if (!q->tip)
	{
		firstraverse(q->next->back, r, sum);
		firstraverse(q->next->next->back, r, sum);
	}
	q->processed = true;
	if (q->tip)
		return;
	kitsch_node * u = q->next->back;
	kitsch_node * v = q->next->next->back;
	unsigned i = u->index;
	unsigned j = v->index;
	unsigned k = q->index;
	if (u->w[j] + v->w[i] <= 0.0)
		q->t = 0.0;
	else
		q->t = (u->w[j] * u->d[j] +
				v->w[i] * v->d[i]) /
			   (2.0 * (u->w[j] + v->w[i]));
	q->weight = u->weight + v->weight + u->w[j] + v->w[i];
//	if (!negallowed)
	assert(q);
   	scrunch(*q);
	secondtraverse(r,q,u,v,i,j,k,sum);
}  /* firstraverse */




//________________________________________________________________________________________

//	sumtraverse

//________________________________________________________________________________________
void sumtraverse(kitsch_node *q, double *sum)
{
	debug_func(3);
	
	/* traverse to finish computation of sum of squares */
	assert(!q->tip);
	if (!q->next->back->tip)
		sumtraverse(q->next->back, sum);
	if (!q->next->next->back->tip)
		sumtraverse(q->next->next->back, sum);
	kitsch_node * u = q->next->back;
	kitsch_node * v = q->next->next->back;
	unsigned i = u->index;
	unsigned j = v->index;
	double TEMP = u->d[j] - 2.0 * q->t;
	double TEMP1 = v->d[i] - 2.0 * q->t;
	(*sum) += u->w[j] * (TEMP * TEMP) + v->w[i] * (TEMP1 * TEMP1);
	//cout << "t sum=" << *sum << "\n";
}  /* sumtraverse */




//________________________________________________________________________________________

//	evaluate

//________________________________________________________________________________________
void evaluate(kitsch_node& r)
{
	/* fill in times and evaluate sum of squares for tree */
	/* Local variables for evaluate: */
	//cerr << "evaluate\n";
	for (unsigned i = 0; i < cnt_nodes; ++i)
		treenode[i]->processed = treenode[i]->tip;
	
	double sum = 0.0;
	firstraverse(&r, &r, &sum);
	if (!r.tip)
		sumtraverse(&r, &sum);
	cnt_trees_examined++;
//	if (replicates && (lower || upper))
//		sum /= 2;
	like = -sum;
	//cout << "e sum=" << sum << "\n";
	//cerr << "evaluate finish\n";
}  /* evaluate */





//________________________________________________________________________________________

//	tryadd

//________________________________________________________________________________________
void tryadd(kitsch_node *p,  kitsch_node** item,  kitsch_node** nufork)
{
	/* temporarily adds one fork and one tip to the tree.
	   if the location where they are added yields greater
	   "likelihood" than other locations tested up to that
	   time, then keeps that location as there */
	//cerr << "tryadd\n";
	add(p, *item, *nufork);
	assert(root);
	evaluate(*root);
	if (like > bestyet)
	{
		bestyet = like;
		there = p;
	}
	re_move(*item, *nufork);
	//cerr << "tryadd finish\n";
}  /* tryadd */




//________________________________________________________________________________________

//	addpreorder

//________________________________________________________________________________________
void addpreorder(kitsch_node *p, kitsch_node* item, kitsch_node* nufork)
{
	/* traverses a binary tree, calling PROCEDURE tryadd
	   at a kitsch_node before calling tryadd at its descendants */
/* Local variables for addpreorder: */
	//cerr << "addpreorder\n";
	if (p == 0)
		return;
	tryadd(p, &item,&nufork);
	if (!p->tip)
	{
		addpreorder(p->next->back, item, nufork);
		addpreorder(p->next->next->back, item, nufork);
	}
	//cerr << "addpreorder finished\n";
}  /* addpreorder */




//________________________________________________________________________________________

//	tryrearr

//________________________________________________________________________________________
void tryrearr(kitsch_node *p,  kitsch_node** r, bool *success)
{
	//cerr << "tryrearr\n";
	/* evaluates one rearrangement of the tree.
	   if the new tree has greater "likelihood" than the old
	   one sets success := TRUE and keeps the new tree.
	   otherwise, restores the old tree */
	kitsch_node *frombelow, *whereto, *forknode;

	if (p->back == 0)
		return;
	forknode = treenode[p->back->index];
	if (forknode->back == 0)
		return;
	double oldlike = like;
	if (p->back->next->next == forknode)
		frombelow = forknode->next->next->back;
	else
		frombelow = forknode->next->back;
	whereto = forknode->back;
	re_move(p, forknode);
	add(whereto, p, forknode);
	if ((*r)->back != 0)
		*r = treenode[(*r)->back->index];
	assert(*r);
	evaluate(**r);
	if (like > oldlike)
	{
		bestyet = like;
		*success = true;
		return;
	}
	re_move(p, forknode);
	add(frombelow, p, forknode);
	if ((*r)->back != 0)
		*r = treenode[(*r)->back->index];
	like = oldlike;
	//cerr << "tryrearr finish\n";
}  /* tryrearr */




//________________________________________________________________________________________

//	repreorder

//________________________________________________________________________________________
void repreorder(kitsch_node *p,  kitsch_node** r, bool *success)
{
	//cerr << "repreorder\n";
	/* traverses a binary tree, calling PROCEDURE tryrearr
	   at a kitsch_node before calling tryrearr at its descendants */
	if (p == 0)
		return;
	tryrearr(p,r,success);
	if (!p->tip)
	{
		repreorder(p->next->back,r,success);
		repreorder(p->next->next->back,r,success);
	}
	//cerr << "repreorder finish\n";
}  /* repreorder */




//________________________________________________________________________________________

//	rearrange

//________________________________________________________________________________________
void rearrange(kitsch_node ** r_)
{
	/* traverses the tree (preorder), finding any local
	   rearrangement which decreases the number of steps.
	   if traversal succeeds in increasing the tree's
	   "likelihood", PROCEDURE rearrange runs traversal again */
/* Local variables for rearrange: */
	//cerr << "rearrange\n";
	kitsch_node **r = r_;
	bool success;
	do
	{
		success = false;
		repreorder(*r,r,&success);
	} while (success);
	//cerr << "rearrange finish\n";
}  /* rearrange */







//________________________________________________________________________________________

//	dtraverse

//________________________________________________________________________________________
/*void dtraverse(kitsch_node *q)
{
	// print table of lengths etc.
	unsigned i;

	if (!q->tip)
		dtraverse(q->next->back);
	if (q->back != 0)
	{
		fprintf(outfile, "%4ld  ", q->back->index - cnt_seq);
		if (q->index <= cnt_seq)
		{
			for (i = 0; i < namelength; i++)
				putc(names[q->index][i], outfile);
		}
		else
			fprintf(outfile, "%4ld      ", q->index - cnt_seq);
		fprintf(outfile, "%13.5f", treenode[q->back->index]->t - q->t);
		fprintf(outfile, "%15.5f\n", root->t - q->t);
	}
	if (!q->tip)
		dtraverse(q->next->next->back);
}  dtraverse */




//________________________________________________________________________________________

//	describe

//________________________________________________________________________________________
void describe(double& sum_sqrs, double& avg_percent_stddev)
{
	/* prints table of lengths, times, sum of squares, etc. */
	

//	fprintf(outfile, "\nSum of squares = %10.3f\n\n", -like);
	sum_sqrs = -like;
//	if (fabs(power - 2) < 0.01)
	{
		double totalnum = 0.0;
		for (unsigned i = 0; i < cnt_seq; ++i)
		{
			for (unsigned j = 0; j < cnt_seq; ++j)
				if (i != j && treenode[i]->d[j] > 0.0)
					totalnum += treenode[i]->w[j] * treenode[i]->d[j] * treenode[i]->d[j];
		}
		totalnum -= 2;
//		if (replicates && (lower || upper))
//			totalnum /= 2;
//		fprintf(outfile, "Average percent standard deviation =");
//		fprintf(outfile, "%10.5f\n\n", 100 * sqrt(-(like / totalnum)));
		avg_percent_stddev = 100 * sqrt(-(like / totalnum));
	}
//	fprintf(outfile, "cnt_trees_examined %4ld trees\n\n", cnt_trees_examined);
//	fprintf(outfile, "From    To           Length          Time\n");
//	fprintf(outfile, "----    --           ------          ----\n\n");
//	dtraverse(root);
//	putc('\n', outfile);
}  /* describe */





//________________________________________________________________________________________

//	copy_kitsch_tree

//________________________________________________________________________________________
void copy_kitsch_tree(vector<kitsch_node*> from, vector<kitsch_node*> to)
{
	/* make a copy of a tree */
	//cerr << "copy_kitsch tree begin\n";
	for (unsigned i = 0; i < cnt_seq; ++i)
	{
		*to[i] = *from[i];
		if (from[i]->back != 0)
		{
			if (from[i]->back == from[from[i]->back->index])
				to[i]->back = to[from[i]->back->index];
			else if (from[i]->back == from[from[i]->back->index]->next)
				to[i]->back = to[from[i]->back->index]->next;
			else
				to[i]->back = to[from[i]->back->index]->next->next;
		}
		else to[i]->back	= 0;
	}
	//cerr << "copy_kitsch tree 1\n";
	for (unsigned i = cnt_seq; i < cnt_nodes; ++i)
	{
		kitsch_node *p = from[i];
		kitsch_node *q = to[i];
		for (unsigned j = 1; j <= 3; ++j)
		{
			*q = *p;
			if (p->back)
			{
				if (p->back == from[p->back->index])
					q->back = to[p->back->index];
				else if (p->back == from[p->back->index]->next)
					q->back = to[p->back->index]->next;
				else
					q->back	= to[p->back->index]->next->next;
			}
			else
				q->back	= 0;
			p = p->next;
			q = q->next;
		}
		//cerr << "!";
	}
	//cerr << "copy_kitsch tree end\n";
}




//________________________________________________________________________________________

//	maketree

//________________________________________________________________________________________
void maketree(std::vector<boost::shared_array<double> >& matrix, double& best_likelihood)
{
	/* constructs a binary tree from the pointers in treenode.
	   adds each kitsch_node at location which yields highest "likelihood"
	   then rearranges the tree for greatest "likelihood" */
	
	for (unsigned i = 0; i < cnt_seq; ++i)
		enterorder[i] = i;
	for (unsigned i = 0; i < cnt_seq; ++i)
		std::swap(enterorder[i], enterorder[rand() % cnt_seq]);
	
	root = treenode[enterorder[0]];
	add(treenode[enterorder[0]], treenode[enterorder[1]], treenode[cnt_seq]);
	//cerr << ".";
	bool lastrearr = false;
	for (unsigned i = 3; i <= cnt_seq; ++i)
	{
		bestyet = -10000000.0;
		kitsch_node* item = treenode[enterorder[i - 1]];
		kitsch_node * nufork = treenode[cnt_seq + i - 2];
		addpreorder(root, item, nufork);
		add(there, item, nufork);
		like = bestyet;
		rearrange(&root);
		assert(root);
		evaluate(*root);
		cnt_trees_examined--;
		lastrearr = (i == cnt_seq);
		//cerr << "4>" << i << "/" << cnt_seq << "<4";
		if (lastrearr)
		{
			//cerr << "5";
			double bestlike = bestyet;
			double gotlike;
			do
			{
				//cerr << "-";
				gotlike = bestlike;
				for (unsigned j = 0; j < cnt_nodes; ++j)
				{
					//cerr << ">";
					there = root;
					bestyet = -32000.0;
					item = treenode[j];
					if (item != root)
					{
					
						re_move(item, nufork);
						there = root;
						addpreorder(root, item, nufork);
						add(there, item, nufork);
					}
				}
			} while (bestlike > gotlike);
			//cerr << "@";
			if (njumble > 1)
			{
				cerr << "bestlike " << bestlike << "\n";
				if (jumb == 1 || (jumb > 1 && bestlike > best_likelihood))
				{
					//cerr << "Z";
					copy_kitsch_tree(treenode, bestree);
					best_root = bestree[root->index];
					best_likelihood = bestlike;
				}
			}
		}
	}
	//cerr << "make tree finish1\n";
	if (njumble == jumb)
	{
		//cerr << "make tree do2\n";
		if (njumble > 1)
		{
			copy_kitsch_tree(bestree, treenode);
			assert(best_root);
			root = treenode[best_root->index];
		}
		//cerr << "make tree evaluate";
		assert(root);
		evaluate(*root);
		//cerr << "make tree finish2";
	}
}  /* maketree */




//________________________________________________________________________________________

//	treeout

//________________________________________________________________________________________
template<typename TREEMAKER> void kitsch_treeout(	TREEMAKER& maker,
													kitsch_node *p, 
													const vector<string>& names)
{
	assert(p);
	if (p->tip)
	{
		maker.add_leaf(names[p->index]);
		// cout << "add_leaf " << p->back<< "\n";
	}
	else
	{
		// cout << "add_internal_node\n";
		maker.add_internal_node();
		kitsch_treeout(maker, p->next->back, names);
		// cout << "add_sibling\n";
		maker.add_sibling();
		kitsch_treeout(maker, p->next->next->back, names);
		// cout << "to_grandparent\n";
		maker.reparent_to_grandparent();
	}
	if (p == root)
	{
		// cout << "finish\n";
		maker.finish();
	}
	else
	{
		// cout << "branch length\n";
		maker.set_branch_len(treenode[p->back->index]->t - p->t);
	}
}
													
string kitsch_tree_to_str(const vector<string>& names)
{
	newhampshire_tree_maker make_newhampshire_tree;
	kitsch_treeout(make_newhampshire_tree, root, names);
	return make_newhampshire_tree.str;
}


//----------------------------------------------------------------------------------------

//	kitsch_make_tree

//----------------------------------------------------------------------------------------
void kitsch_make_tree(	paralogy_tree::t_tree& tree,
						const std::vector<std::string>& names,
						std::vector<boost::shared_array<double> >& matrix,
						double& sum_sqrs, double& avg_percent_stddev,
						unsigned& cnt_trees_examined,
						t_progress_indicator& dots)
{
	cnt_seq = matrix.size();
	cnt_nodes = cnt_seq * 2 - 1;
	//cout << "init\n";
	doinit();
	enterorder.insert(enterorder.begin(), cnt_seq, 0);
	assigndata(matrix, cnt_seq, cnt_nodes);
	cnt_trees_examined = 0;
	double best_likelihood = 0.0;
	for (jumb = 1; jumb <= njumble; jumb++)
		maketree(matrix, best_likelihood);
	
	// use best_root tree
	if (njumble > 1)
	{
		copy_kitsch_tree(bestree, treenode);
		root = treenode[best_root->index];
	}
	evaluate(*root);
	
	//cout << "describe\n";
	describe(sum_sqrs, avg_percent_stddev);


	//cout << "make_paralogy_tree\n";
	paralogy_tree_maker make_paralogy_tree(tree);
	kitsch_treeout(make_paralogy_tree, root, names);
	//cout << "clean up memory\n";
	enterorder.clear();
	shared_mem.clear();
}

//________________________________________________________________________________________

//		kitsch_make_tree

//			wrapper

//________________________________________________________________________________________
void kitsch_make_tree(	paralogy_tree::t_tree& tree,
						const std::vector<std::string>& names,
						std::vector<boost::shared_array<double> >& matrix,
						double& sum_sqrs, double& avg_percent_stddev,
						unsigned& cnt_trees_examined)
{
	ostringstream ignore_strm;
	t_progress_indicator dots(ignore_strm, 10000);
	kitsch_make_tree(tree, names, matrix, sum_sqrs, avg_percent_stddev, 
						cnt_trees_examined, dots);
}






}

